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On  Curve  Interpolation  in  IRd 


Jernej  Kozak  and  Emil  Zagar 


Abstract.  In  this  paper  the  interpolation  by  G 2 continuous  spline  curves 
of  degree  n in  lRd  is  studied.  There  are  r interior  and  two  boundary  data 
points  interpolated  on  each  segment  of  the  spline  curve.  The  general  form 
of  the  spline  curve,  as  well  as  the  defining  system  of  nonlinear  equations 
are  derived.  The  asymptotic  existence  of  the  solution,  and  the  approxi- 
mation order  are  studied  for  the  polynomial  case  only.  It  is  shown  that 
the  optimal  approximation  order  is  achieved,  and  asymptotic  existence  is 
established  provided  the  relation  r = n — 2 is  satisfied.  These  conclusions 
hold  independently  of  d.  It  is  also  pointed  out  that  the  underlying  analysis 
could  not  be  carried  over  to  the  case  r = n — 1. 


§1.  Introduction 

The  interpolation  problem  considered  is  the  following.  Let  the  points 

Tq,Ti,  . . . ,Tjv  € lRd,  Tj  ^ Tj+i,  all  j,  d > 2,  (1) 

and  the  tangent  directions 

d0,du  (2) 

at  the  boundary  points  be  given.  Find  a G2  continuous  spline  curve  Bn  of 
degree  n which  interpolates  the  prescribed  data. 

The  problem  appeared  first  as  a particular  limit  case  in  [2],  and  was 
further  generalized  in  several  papers,  among  them  in  [3-5,6,9-10].  A general 
approach  to  the  approximation  order  achieved  can  be  found  in  [8] . 

Here,  the  general  setup  is  tackled.  The  interpolating  spline  curve  in  the 
Lagrange  form  is  established  and  the  defining  system  of  nonlinear  equations 
is  derived  in  general.  However,  the  asymptotic  existence  of  the  solution  (i.e. 
the  existence  of  the  solution  when  given  points  are  sampled  densely  enough) 
and  the  approximation  order  turned  out  too  comprehensive  to  be  studied  here 
in  a general  framework.  The  positive  conclusions  for  the  single  segment  case 
when  r = n — 2 are  established.  It  is  possible  to  extend  these  results  to 

Curve  and  Surface  Fitting:  Saint-Malo  1999  263 

Albert  Cohen,  Christophe  Rabut,  and  Larry  L.  Schumaker  (eds.),  pp.  263-272. 

Copyright  ©2000  by  Vanderbilt  University  Press,  Nashville,  TN. 

ISBN  0-8265-1357-3. 

All  rights  of  reproduction  in  any  form  reserved. 


264 


J.  Kozak  and  E.  Zagar 


the  m-segment  spline  curve,  but  the  proofs  are  not  short,  and  will  appear 
elsewhere.  On  the  contrary,  as  one  could  guess  from  [8],  the  case  r — n — 1 is 
not  encouraging. 

Why  would  one  use  the  G2-continuous  splines  as  interpolating  curves? 
Quite  clearly,  the  derivative  continuity  at  the  breakpoints  becomes  in  this 
way  independent  of  the  local  parametrisation.  Also,  these  curves  could  be 
seen  as  a generalization  of  the  odd  order  spline  function  interpolation  at  knots, 
applied  so  successfully  in  many  cases.  The  order  of  G-continuity  2 is  pinned 
down  by  the  human  eye,  sometimes  quite  important  in  CAGD:  it  can  detect 
the  continuity,  the  continuity  of  the  tangent  direction  and  the  curvature,  but 
hardly  higher  order  geometric  quantities. 

Throughout  the  paper  bold  faced  letters  will  stand  for  vectors,  and  or- 
dinary ones  for  scalars.  The  dot  product  on  ftd  will  be  denoted  by  • and  its 
implied  norm  by  ||-||.  Derivatives  with  respect  to  the  global  (or  local)  pa- 
rameter will  be  denoted  by'  (or  d/d(),  and  those  with  respect  to  the  natural 
parameter  by  '. 

Now  let  Bn  be  a continuous  spline  curve  of  degree  n with  m segments 
B :=  Bn  : [CoiCm] 
corresponding  to  the  breakpoints 


Co  Cl  < Cm  • 


We  suppose  that  the  pieces  are  locally  parametrized  on  [0, 1]  as 

m = Ce[6-i,6]. 

Suppose  B interpolates  the  data  (1),  and  (2).  If  r interior  and  two  boundary 
points  are  to  be  met  on  each  segment,  then  N = m(r  + 1).  Further,  on  the 
Gth  segment  the  interpolation  conditions  read 

= Ttj  :=  T(/_i)(r+i)+j,  j = 0,  + 1,  C=l,2,  (3) 


where 

0 =:  t(to  < t(  i < • • • < ti,r+ 1 :=  1, 


and  (G,j)j=i  are  the  unknown  parameters  to  be  determined.  Let  x A y >— > 
{xiTjj  — XjVi)i<j  denote  the  2-wedge  product.  The  geometric  continuity  of  B 
requires  the  tangent  direction 


(4) 


1 


BhB 


as  well  as  the  curvature 


(5) 
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d6 


Fig.  1.  An  interpolating  spline  curve  with  three  segments  . 


to  be  continuous  at  the  breakpoints.  Additionally,  at  the  boundary  points  the 
tangent  directions  do  and  have  to  be  interpolated  too,  i.e., 


do  A B( Co)  = B((m)  A d/v  = 0.  (6) 

Fig.  1 gives  an  example  of  such  an  interpolating  spline  curve  for  r = 1,  n = 3, 
and  d = 3.  A brief  look  at  the  conditions  (3)-(6)  reveals  that  the  number  of  in- 
dependent equations  would  be  equal  to  the  number  of  independent  unknowns 
if 


dn  — (d  — 1)  r = 3d  — 2. 


(7) 


As  already  observed  in  [5],  for  fixed  d this  Diophantine  equation  always  has 
an  infinite  number  of  nonnegative  solutions.  The  following  lemma  gives  its 
general  solution. 


Lemma  1.  The  possible  choices  of  pairs  r and  n that  satisfy  (7)  for  fixed  d 
are  given  by 


r = d — 2 + dk,  n = d + (d  — l)k,  fc  = 0, 1, (8) 

Proof:  The  relation  (7)  can  be  rewritten  as 

d(n  — d)  — (d  — l)(r  — d + 2)  = 0. 

Since  d > 2,  the  numbers  d and  d — 1 are  relatively  prime.  So  d must  divide 
r — d + 2,  and  d — 1 must  divide  n — d,  i.e., 

r — d + 2 n — d 

~^—=d^i  = k 

for  an  integer  k.  But  r — d — 2 + dk  > 0 implies  k > | — 1 > — 1,  and  the 
conclusion  follows.  □ 
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§2.  The  Defining  Equations 

Several  approaches  were  used  to  simplify  the  conditions  (3)-(6)  for  particular 
choices  of  d,  n,  and  r.  Here  we  show  that  this  can  be  done  in  general,  which 
will  provide  an  opportunity  to  unify  the  computer  programs.  Let  us  consider 
a single  segment  first.  In  this  case,  the  data  to  be  interpolated  are  the  points 
T0,T1,...,Tr+1,  Ti?  Tj+ 1,  as  well  as  tangent  directions  do,  dr+j  at  the 
boundary  points.  Suppose  r and  n are  given  by  (8).  Consider  the  case  n = r+2 
first,  i.e.,  k = 0.  The  interpolating  polynomial  curve  can  be  written  explicitly 
in  Lagrange  form  as 

r+l 

B b uj  + ^ ' T j£j 

1=0 

with 

:=  rp  - :=  (t ~ tj)l'(tj)’  (9) 

tj  :=  tlfj , and  the  values  (tj)j=1,  to  be  determined.  Here  b £ Rrf  denotes  the 
unknown  leading  coefficient  vector.  If  k > 1,  one  has 

r + 2 = d(k  + 1)  > d(k  + l)  - k = n, 

and  B is  of  degree  at  most  r + l,  i.e., 

r + l 

B = J2TiCr 

i=o 

In  particular,  this  imposes  additional  conditions 

r+l 

degree  ^ Tj  Cj  < n (10) 

l=o 

for  k > 1.  An  easy  way  to  meet  the  tangent  direction  conditions  (6)  is  to 
introduce  two  additional  (strictly  positive)  real  unknowns,  Qo  and  ar+i,  and 
require 

B(to)  = aodo,  B(tr+i)  = ar+idr+i.  (11) 

Let 

t_i  = To  :=  <o>  ij'  tj  j j — 1)  2, . . . , r,  rr+2  rr_)_i.  tr+1  * (12) 

Since  B is  a polynomial  of  degree  < n,  the  divided  difference,  based  upon 

n + 2 — r + A — k 

points  maps  it  to  zero.  So  the  conditions  (11)  and  (10)  can  be  written  in  a 
compact  form  as 

[rj-i,  Tj, . . . ,rJ+r+2-fc]S  = 0,  j = 0, 1,  — , fc,  (13) 
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which  is  a system  of  d(k  + 1)  nonlinear  equations  for  r + 2 = d(k  + 1)  scalar 
unknowns 

0^0 1 ^1 ) ^2 ) • • ■ j j ^r+l  • (1^) 

In  the  case  n = r + 2,  one  has  to  determine  additionally  the  coefficient  vector 


b,  for  example  as 

b = [t0 , to,  1 1 j . . . , tr , tr+i\B  — [to,  t\ , . . . , tT , tr_|_i , tr+i]B . (1 5) 

Now,  for  an  m-segment  spline  curve,  the  directions  d( , i — 1, 2, . . . , m — 1,  are 
unknown,  as  well  as 

^t,0,  t^l,  t^( 2)  • • * i t^r,  0^}r+ 1 , ^ 1,2,...,  TCI.  (10) 

But  one  can  still  write  the  interpolation  conditions  on  the  f-th  segment  as 

1 > * * • i Tgj+r+2— k\B  “0,  j 0, 1, . . . , k,  (If) 


where  tij  are  defined  as  in  (12),  but  this  time  for  the  composite  case.  In 
addition,  the  missing  ( d — l)(m  — 1)  equations  are  supplied  by  the  continuity 
conditions  of  the  curvature  (5). 

§3.  Asymptotic  Existence  and  Approximation  Order 

The  system  of  equations  based  on  (17)  and  continuity  of  curvature  (5)  is  non- 
linear, and  one  of  the  approaches  to  study  it  is  to  assume  that  the  data  (1)  and 
(2)  are  based  upon  a smooth  underlying  regular  parametric  curve  / : I — > lllrf, 
parametrized  by  the  arclength  s.  The  local  expansion  of  the  curve  /,  and  the 
data  T(j  (sampled  densely  enough),  give  rise  to  an  asymptotic  analysis  of  the 
nonlinear  system.  The  simplest  way  to  obtain  the  local  expansion  is  to  use  the 
Frenet  frame  as  the  local  coordinate  system,  and  the  Frenet-Serret  formulae 
to  obtain  this  expansion.  Let  (e,(s))f=1  denote  the  Frenet  frame,  with 

/'  = «!•  (18) 


The  Frenet-Serret  formulae  read 
ei(s)  = Ki(s)e2(s), 

e'(s)  = -/ti_i(s)ej_i(s)  + Ki(s)ei+1(s),  i = 2,  3, . . . , d - 1,  (19) 

ed(«)  = -Kd-i(s)ed_i(s), 


where  Kj  are  first  d — 1 principal  curvatures  of  /,  expanded  as 

. . 1 1 2 
Kj(s)  = KiO  + + yKi2S  +■■■ 


(20) 


Since  / is  a regular  curve,  k,,o  > const  > 0,  i = 1,2, ...,d  — 2.  We  will 
additionally  assume  that  k^o  > const  > 0.  Beginning  with  (18),  the  higher 
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derivatives  of  / can  be  computed  by  (19)  and  (20).  This  produces  the  required 
expansion 

/W  = /(0)  + /'(0)s  + i/"(0)S2  + -.. 

= /(0)  + (s-  jU?,0s3  + -'-)ei(0)  (21) 

+ (x«l,0S2  + X«l,lS3  + • ■ -)e2(0)  + (^K1,0'C2,0S3  + • • -)e3(0)  + • • • 

2 b b 

Let  us  now  consider  the  single  segment  case  of  the  interpolation  problem 
with  data  based  on  a smooth  / : [0,  h]  — * ]Rd, 

do  = f'{mh),  Tj  = f(r)jh),  j = 0,1,..., r + 1,  dr+l  = f'{r)r+1h), 

with  points  separated  independently  of  h,  i.e., 

0 =:  7?0  < jji  < • • • < < r)r+\  :=  1. 

Since  translation  and  rotation  do  not  influence  the  asymptotic  analysis,  we 
may  assume  /( 0)  = 0,  and 

*(0)  = (Mi=i.  * = (22) 

Then,  with  the  help  of  (21),  one  obtains 

f(Vjh)=(^r,}hi,f[Kq,0  + O(hi+1))  , (23) 

V ' «=o  ) <=i 

and  a similar  expression  for  Since  the  divided  difference  is  a linear 

functional,  we  can  normalize  the  system  (13)  by  multiplying  the  data  values 
by  D-1,  with 

(l  . <_1 

D :=  diag  I ^ h’  kQi0 

Let  f(s)  :=  (sl)f=1  denote  the  leading  part  of  the  normalized  /.  Then 

[to,to,tu  ■ ■ ■ , tr,tr+i,tr+i]D  lB  — [to  ito,ti, . . . ,fr,tr_|.i,fr-|-i].Z?  + O(h), 

(24) 

and  B is  a polynomial  of  degree  < n = r + 2 that  satisfies  the  interpolation 
conditions 

B (tj)  = OLjf  (r}j),  i = 0,r  + l, 

B(tj)  — J(Vj)>  J =0,l,...,r,r  + 1, 

~ ao  ~ ar+i 
a0:-T,  ar+1:=—. 


where 
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Note  that  all  the  components  of  / are  polynomials  of  degree  < d = r + 2.  This 
implies  that 

[»Zo,  »7o,  >71,  • • • , Vr,Vr+i , Vr+i]f  = 0,  (25) 

and  the  solution  of  (24)  in  the  limit  h —►  0 now  reads  as 

t*  :=  (5 5, . . . Xr,K+i)  = (i, m,V2,  — i).  (26) 

To  prove  the  existence  of  the  solution  for  h small  enough,  it  is  sufficient  to 
show  that  the  Jacobian  of  the  system  (24)  is  nonsingular  at  the  limit  (26). 
The  Jacobian  will  be  determined  with  the  help  of  the  following  fact:  if  xj  is 
different  from  all  the  other  points  Xi,  and  if  a function  g is  smooth  enough, 
one  has 


,*>,• 


> xji  • 


n (xi  - xi) 


[. . . ,Xj,Xj, . . ,\g 


g'jxj) 

n — x^ 


(27) 


Consider  now  B = (B  — f)+f.  Since  B— f = 0 at  t*,  all  its  partial  derivatives 
with  respect  to  tj  vanish,  and  this  difference  contributes  to  the  Jacobian  at 
the  limit  point  t*  only  in  the  first  and  last  column,  i.e., 


— [t0,  to,  <1,  • • ■ , fr,  tr+1,  tr+1]  (B  - J)  (%), 


d 1 

w [to,  to,  ti, . . . , tr,  tr_|_i,  tr+i]  \B  f ) | = v~77  V 

ocer+ 1 It * (Vt+ l-lloKWr+l) 


where  to  is  given  by  (9),  and 


/ (»?r+l)i 
(28) 


U)  U)\ 


The  polynomial  curve  / does  not  depend  on  oto,  av+i,  and  from  (27)  and  (25) 
one  obtains  the  columns  2, 3, . . . , r + 1 with  j = 1, 2, . . . , r as 


[to,  to,  tl,  . . . , tr,  tr_|_i,  tr+l] 


/ It* 


{Vj  - Vo ){Vj  - Vr+l)dj'{tj) 


f (Vj)- 


It  is  now  straightforward  to  see  that  the  Jacobian  at  t*  is  the  Vandermonde 
matrix  V(?7o,*?i>  • • • , Vr+i),  multiplied  by  D i :=  diag(i)f=1  from  the  left,  and 

by 


I>2  :=  diag 


Oj'ivo)  ’ »7l(l  - Vl  V(»?l)  ’ ' " ’ Vr(l  - VrWiVr)  ’ ^{Vr+l) 


from  the  right.  This  prepares  the  proof  of  the  following  theorem. 
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Theorem  2.  The  system  (13)  has  a unique  solution  for  h small  enough.  The 
approximation  order  of  the  resulting  interpolating  polynomial  curve  Bn  is 
optimal,  i.e.,  r + 4 = n + 2. 

Proof:  Since  the  matrices  P(r/o,»?i, . . . , J?r+i),  D\  and  D2  are  nonsingular, 
the  Jacobian  at  the  limit  point  t*  is  nonsingular,  too,  and  the  existence  of  a 
unique  solution  for  h small  enough  is  established.  Furthermore,  the  unknown 
parameters  are  of  the  form 

a0  = ar+ 1 = h + 0(h2),  tj  = rjj  + 0(h),  j = 1, 2, . . . , r.  (29) 

Since  there  are  r + 2 points,  as  well  as  two  directions  interpolated,  the  optimal 
approximation  order  is  quite  clearly  < r + 4.  The  proof  will  now  follow  the 
approach  applied  in  [2],  and  extended  in  [5].  It  is  based  on  a reparametrisation 
that  transforms  the  direction  interpolation  to  the  derivative  interpolation,  and 
gives  an  estimate  of  the  parametric  approximation  order  as  defined  in  [7]. 
Recall  (22),  and  the  fact  that  interpolation  conditions  are  satisfied.  By  [2] 
and  [5],  it  is  now  enough  to  confirm  that  all  the  components  of  / and  B can 
be  reparametrized  by  the  ordinate  of  the  first  component  of  both  curves.  As 
to  /,  for  h small  enough  this  fact  is  obvious.  The  first  component  behaves 
by  (21)  as  s + 0(s3),  and  the  others  at  least  as  0(s2).  To  establish  the  same 
conclusion  for  B,  it  is  enough  to  show  that 

B - ch(Su)f=1  + 0(h2),  c 7^  0.  (30) 

Further,  the  optimal  approximation  order  proof  depends  on  the  additional 
relations 

B^=0(hq),  q = 2, 3,...,  r + 2.  (31) 

The  result  required  then  follows  from  the  standard  error  estimate  of  interpo- 
lation, and  the  fact  that  the  (r  +4)-th  derivative  of  B with  respect  to  the  new 
parameter  is  bounded  independently  of  h.  Let  us  verify  the  relations  (30)  and 
(31).  Recall  first 

^EW),  <7  = 0,l,...,r  + l,  tr+2=w(t)  + ^t;+%(<).  (32) 

3=  0 3=0 

The  divided  difference  [to,  to,  t\, . . . , tr,  tr+i]  maps  polynomials  of  degree  < 
r + l = d— lto  zero,  and  depends  continuously  on  its  arguments  if  applied  to 
a smooth  function.  Thus  b by  (15)  and  (23)  near  the  limit  point  t*  behaves 
like 

b = ( o(hd ),  o(hd), . . . , o(hd),  Xd  hd  + 0(hd+1))T , 
where  Xi  = 11^=0  Kq, 0 > 0-  On  !he  other  hand,  (29)  and  (32)  imply  that 

^TiCj(t)  = ;£f(Vjh)Cj(t) 

j= 0 3=0 

= (xi  ht,...,  Xd-1  h "-1  td~\  Xd  hd  (td  - w(t)))T  (1  + 0(h)). 
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But 

r+l 

B(,)(t)  = b w(*)W  + J2  q = 1, 2, ....  r + 2, 

j=o 

and  (31)  follows.  The  proof  is  complete.  □ 

There  is  no  hope  that  this  approach  could  be  used  for  all  k.  In  fact,  it 
fails  already  for  k — 1,  as  we  will  show  now.  By  (13),  the  equation  (24)  is 
replaced  by 

• • • >tr>tr+l]D  1B  = [to,  tO)£l>  • • ■ )tr,  U+i]B  + 0(h), 

[fo,ti, . . . ,tr,tr+i,tr+i\D~1B  = [t0,ti, . . . ,tT,tT+i,tr+i]B  + 0(h). 

Further,  as  in  the  proof  of  Theorem  2,  the  first  column  of  the  Jacobian  is 
determined  from 

— [t0,to,ti,...,tr,tr+i](B-7)j^  = ~7^y 7 M, 

3 

^0;  tl,  • • • , tr , tr+li  tr+l]  {B  — f)  | 0, 

the  last  column  from 

■ j ^ ^ 

[to,ta,ti,...,tr,tr+i](B-/)|  =0, 

<JOtr+ 1 It* 

[t0,  tl,  • • • , tr,  tr+l,  tr+l]  ( B - f)  (ljr+l), 

and  the  other  columns  from 

(^-[t°,t0,t1,...,tr,tr+1])  7|f  = 

(st;ltO,tl,-"’*r,tr+1,tr+10  = ~(Vj  -Tlr+J&fa)* 

After  normalizing  the  Jacobian  from  the  left  by  Df1,  and  by  D^1  from  the 
right  one  obtains  the  matrix  A :=  (ay)ij=i  with 

^i,i  = ^i,i,  t — 1,2,...,  2d, 

& i,2d  ” 0,  ^i-\-d,2d  — 1,  t — 1,  2,  . . . , d, 

and 

®-i,j  — Vj_ i Vj—ii  ®i+d,j  Vj—h  t 1, 2, . . . , d,  j 2,3,...,  2d  1. 

A simple  rank  preserving  transformation 

^ — i,j,  t — 2d,  2d  1, . . . , d T 1,  ) 1, 2, . . . , 2d, 

transforms  A to  a matrix  with  row  i equal  to  row  i + d for  i = 2, 3, . . . , d.  It 
is  now  easy  to  see  that  the  rank  of  the  matrix  A is  d + 1,  and  consequently 
dim  ker  A = d — 1.  Thus,  since  the  Jacobian  is  singular,  some  other  approach 
such  as  [1],  pp.  154-155,  should  be  applied  to  carry  out  the  asymptotic  anal- 
ysis. 
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